#!/usr/bin/env bash

###check for parameters' number
if [ $# -lt 4 ];then
	echo -e "error : date_of_freqfiles liftover.sh mergefile databasefile file1 file2 ... \n file1 = liftover_wellconverted_dataname1.txt, file2 = liftover_wellconverted_dataname2.txt etc..." 1>&2
	exit 1
fi

### check for frequencies files date
if [ -d $1 ];then
	echo "The first parameter (date_of_freqfiles) must be a date number(e.g. frequencies_dataname_1232011.out, date_of_freqfiles : 	1232011)"
	exit 2
else
	if [ $(($1+0)) -eq 0 ];then # eq 0 if it isn't a number
		echo "The first parameter (date_of_freqfiles) must be a date number(e.g. frequencies_dataname_1232011.out, date_of_freqfiles : 	1232011)"
	exit 2
	else
		datefreqfile=$1
	fi
fi
shift

datefile=`date +%m%e%Y` # date mmddyyyy

### mergefile
echo -e "\n************* CHOOSE OR LOAD MERGE FILE *************"
find merge_[0-9]*.txt > /dev/null
if [ $? -ne 0 ];then
	echo "load $1"
	cat $1 | tr '	' ' ' | sort -t' ' -k1,1 > "merge_"$datefile".txt"
	mergefile="merge_"$datefile".txt"
else
	echo -e "\nmerge file(s) in your directory : "
	find merge_[0-9]*.txt 
	echo -e "\nenter merge filename to use or press \"c\" to load merge parameter's file"
	read choice1
	if [ "$choice1" == "c" ];then
		echo "load $1"
		cat $1 | tr '	' ' ' | sort -t' ' -k1,1 > "merge_"$datefile".txt" 
		mergefile="merge_"$datefile".txt"
	else
		mergefile=$choice1
	fi
fi
shift

### databasefile
echo -e "\n*********** CHOOSE OR LOAD DATABASE FILE ************"
find db_[0-9]*.txt > /dev/null
if [ $? -ne 0 ];then
	echo "load $1"
	cat $1 |  awk -F" " '{print $1"-"$3,$1,$2,$4,$5}' | sort -t' ' -k1,1 > "db_"$datefile".txt" #sort 1000G database file by chromosome and position
	db="db_"$datefile".txt"
else
	echo -e "\ndatabase(s) file(s) in your directory : "
	find db_[0-9]*.txt 
	echo -e "\nenter database filename to use or press \"c\" to load database's parameter's file"
	read choice2
	if [ "$choice2" == "c" ];then
		echo "load $1"
		cat $1 |  awk -F" " '{print $1"-"$3,$1,$2,$4,$5}' | sort -t' ' -k1,1 > "db_"$datefile".txt" #sort 1000G database file by chromosome and position
		db="db_"$datefile".txt"
	else
		db=$choice2
	fi
fi
shift

### create output directory if it doesn't exist yet
if [ ! -d "./output_liftoversh" ]; then
	mkdir output_liftoversh
fi

echo -e "\n************** START FILES' TREATMENT ***************\n"
for f in $@;do

	echo "FILE : $f " 
	### datafiles wellconverted "rs" SNPs
	name=`basename $f .txt` # remove ".txt" suffix
	cat $f | sed 's/^chr//' | awk -F"	" '{print $1"-"$2,$4}' | sort -t' ' -k1,1 > "sort_"$name"_"$datefile".txt" # sort by chromosome and position
	echo "chr-position	snpIDdata	chr	snpIDdb		a0DB	a1DB" > "output_liftoversh/summary_"$name"_"$datefile".txt" # header
	join -t' ' -1 1 -2 1 "sort_"$name"_"$datefile".txt" $db | tr ' ' '	' >> "output_liftoversh/summary_"$name"_"$datefile".txt"  # join files by identical chromosome-position
	echo "see info :  output_liftoversh/summary_"$name"_"$datefile".txt"
	
	# select lines with rsID in 1000genomes database 
	echo "chr-position	snpIDdata	chr	snpIDdb		a0DB	a1DB" > "output_liftoversh/rs_"$name"_"$datefile".txt"  # header 
	cat "output_liftoversh/summary_"$name"_"$datefile".txt" | egrep "	[0-9]+	rs" >> "output_liftoversh/rs_"$name"_"$datefile".txt" 

	echo "see info : output_liftoversh/rs_"$name"_"$datefile".txt"

	### datafiles wellconverted "chr-pos" SNPs and frequencies' files
	if [ ! -f "frequencies"${name#liftover_wellconverted}".bim_"$datefreqfile".out" ];then
		echo "for frequencies_dataname.bim_dateoffreqfile.out liftover file must be : liftover_wellconverted_dataname.txt and liftover_failed_dataname.txt"
	exit 3
	fi
	cat "frequencies"${name#liftover_wellconverted}".bim_"$datefreqfile".out" | sort -t'	' -k2,2 | awk -F"	" '{print $2,$6,$7,$8,$4,$5}' | tr '	' ' ' > tmpfreq.txt
	# select lines with chrposID in 1000genomes database 
	cat "output_liftoversh/summary_"$name"_"$datefile".txt" | egrep -v "	[0-9]+	rs" | sort -t'	' -k2,2 | tr '	' ' ' > "tmpchrpos_"$name"_"$datefile".txt" 
	echo "chr-position	snpIDdata	chr	snpIDdb		a0DB	a1DB	freqDB_a0 freqdata_a0	deltafreq	a0data	a1data" > "output_liftoversh/chrpos_"$name"_"$datefile".txt" #header of chrpos file
	join -1 2  -2 1  "tmpchrpos_"$name"_"$datefile".txt" tmpfreq.txt | awk -F" " '{print $2,$1,$3,$4,$5,$6,$7,$8,$9,$10,$11}' | tr ' ' '	' >> "output_liftoversh/chrpos_"$name"_"$datefile".txt"  

	echo "see info : output_liftoversh/chrpos_"$name"_"$datefile".txt"
	
	rm "sort_"$name"_"$datefile".txt" 

	### datafiles failed
	endfailedfile=`basename ${f#liftover_wellconverted} .txt`
	failedfile="liftover_failed"$endfailedfile
	cat $failedfile".txt" | sed '/#Deleted in new/d' | cut -d"	" -f4 | uniq > "output_liftoversh/exclude_"$failedfile"_"$datefile".txt"
	echo "see list of SNPs to del : output_liftoversh/exclude_"$failedfile"_"$datefile".txt"

	### check RsMerge
	# delete "rs" prefix 
	cat "output_liftoversh/rs_"$name"_"$datefile".txt" | sed ':z;s/rs//;tz' > temporary.txt
	# search rsID 1000G in RsMerge
	cut -d"	" -f4 temporary.txt | sort > listotest.txt
	join -1 1 -2 1 listotest.txt $mergefile > "jointmp2_"$name"_"$datefile".txt"
	cat "jointmp2_"$name"_"$datefile".txt" | sort -t' ' -k2,2 > "jointmp3_"$name"_"$datefile".txt"
	cut -d"	" -f2 temporary.txt | sort > listotest2.txt
	join -1 1 -2 2 listotest2.txt "jointmp3_"$name"_"$datefile".txt" | tr ' ' '	' > "tmpkeep_"$name"_"$datefile".txt"
	cat "tmpkeep_"$name"_"$datefile".txt" | cut -d"	" -f2 | awk '{print "rs"$1}' > "output_liftoversh/keep_"$name"_"$datefile".txt"
	echo -e "see list of SNPs to keep : output_liftoversh/keep_"$name"_"$datefile".txt \n" #-e allows \n interpretation

done

echo -e "\n*********************** ALL DONE *************************"

# delete temporary files
rm tmpfreq.txt
rm tmpchrpos_*
rm temporary.txt
rm listotest.txt
rm listotest2.txt
rm jointmp*
rm tmpkeep*

# command example : ./liftover3.sh 7212011 ../../db/RsMergeArch.bcp ../../db/NEW_1000/EUR.chrfulldata.db liftover_wellconverted_*
